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We show that radiation from complex and inherently random but correlated wave sources can be 
modelled efficiently by using an approach based on the Wigner distribution function. Our method 
exploits the connection between correlation functions and the Wigner function and admits in its sim¬ 
plest approximation a direct representation in terms of the evolution of ray densities in phase space. 

We show that next leading order corrections to the ray-tracing approximation lead to Airy-function 
type phase space propagators. By exploiting the exact Wigner function propagator, inherently 
wave-like effects such as evanescent decay or radiation from more heterogeneous sources as well as 
diffraction and reffections can be included and analysed. We discuss in particular the role of evanes¬ 
cent waves in the near-held of non-paraxial sources and give explicit expressions for the growth 
rate of the correlation length as function of the distance from the source. Furthermore, results for 
the reffection of partially coherent sources from hat mirrors are given. We focus here on electro¬ 
magnetic sources at microwave frequencies and modelling efforts in the context of electromagnetic 
compatibility. 
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I. INTRODUCTION 

Predicting the properties of wave helds in complex en¬ 
vironments is an extremely challenging task of crucial 
importance to a wide variety of technological and en¬ 
gineering applications, such as vibroacoustics [l| or elec¬ 
tromagnetic (EM) wave modelling. In particular, charac¬ 
terising the radiation of EM sources reliably, both in free 
space and within enclosures, is a longstanding research 
issue. In the context of electromagnetic compatibility 
(EMC), digital circuits and large printed circuit boards 
(PCB) embed thousands of electronic devices and metal¬ 
lic tracks and can produce helds reaching dangerous but 
hard-to-predict levels [2|. 

In this paper, we set out an approach for propagat¬ 
ing such complex and statistically characterized wave 
helds exploiting Wigner distribution function (WDE) 
techniques. This approach has its origin in quantum 
mechanics Q, but has more recently found widespread 
attention in optics, see for an overview. The WDE 
formalism offers a direct route to pure ray-tracing ap¬ 
proximations in an operator implementation Q, while 
still capturing in its exact formulation the full wave dy¬ 
namics. The formalism allows one to efficiently treat ra¬ 
diation from complex sources, often having a signihcant 
nondeterministic, statistical character. 

The method introduced below exploits a connection 
between the held-held correlation function (CE) and the 
WDF lillii. Both quantities have been studied inten¬ 
sively in the physics and optics literature. Eor wave 
chaotic systems. Berry’s conjecture postulates a univer¬ 
sal CE equivalent to correlations in Gaussian random 
helds Non-universal corrections can be retrieved 


by linking the CE to the Green function of the system 
[l2l-[l5j. In this paper, we describe how held-held CEs 
can be efficiently propagated using ideas based on ray 
propagation in phase-space. We discuss furthermore non- 
paraxial effects as well as including near held effects due 
to evanescent wave contributions. A systematic expan¬ 
sion of the Wigner function propagator including next-to- 
leading-order effects in the propagating regime leads to 
Airy-function integral kernels containing the ray-tracing 
propagator in the small wavelength limit, akin to the 
treatments in [3l|, [33| . We also show that our WDE rep¬ 
resentation conhrms the validity of the generalised form 
of the Van Cittert-Zernike (VCZT) theorem discussed 
in [TgI. We give a natural extension of this generalised 
VCZT for non-paraxial sources and in the near-held re¬ 
gion where evanescent waves play a prominent role. 

We illustrate these techniques in the context of appli¬ 
cations in EMC and related issues. Here, the system un¬ 
der investigation represents a high-density interconnect 
of integrated electronic circuits. Simulating EM held dis¬ 
tributions in a reliable way is highly topical; in addition, 
the wave correlation function can be measured explic¬ 
itly in this regime 0, thus providing the necessary in¬ 
put information for numerical simulations. The system 
under consideration in this paper consists of a series of 
parallel tracks carrying partially correlated currents, and 
mimics the typically very complex EM sources found on 
PCBs. The method has much wider application, how¬ 
ever. In particular, when combined with fast phase-space 
propagation methods such as the Discrete Flow Mapping 
techniques developed in the context of vibro-acoustics 
Hi, the proposed WDE approach offers an ideal plat¬ 
form for developing a universal high-frequency simulation 
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method. 


II. PHASE-SPACE REPRESENTATION OF 
CLASSICAL FIELDS 

Radiation from simple EM sources such as antennae 
can be characterised deterministically through classical 
electrodynamics methods M- Even though such sources 
are regular and homogeneous, efficiently predicting far- 
held emission from the near-held pattern requires non¬ 
trivial ehort if the sources are extended over many wave¬ 
lengths [T^ . EM sources are becoming increasingly com¬ 
plex, however, and the problem of radiation from dig¬ 
ital circuits or PCBs presents even greater challenges. 
Modelling such sources deterministically is often infeasi¬ 
ble due to the complexity of the structures, whose de¬ 
tails may not even be known in practice. Each compo¬ 
nent of such a complex EM source is typically driven by 
unknown sets of quasi-random voltages, subject to fast 
transients 13. This is due to the presence of a multi¬ 
tude of electronic components whose switching behaviour 
depends on the instantaneous operation mode of the cir¬ 
cuit, and whose excitation signals are intrinsically ran¬ 
dom, or highly sensitive to frequency HH. Consequently, 
the physical investigation of these scenarios challenges 
existing analytical and numerical techniques, and calls 
for more sophisticated modelling tools. 

It is thus natural to use statistics as a language for de¬ 
scribing the radiation from such complex sources. Specif¬ 
ically, we do not attempt to characterise or propagate the 
field itself, which is typically hard to obtain in practice, 
but rather its two-point CE. It has been demonstrated in 
[HI that corresponding measurements are feasible in the 
context of emission from electronic devices and PCBs. 
Here we describe the basic elements needed to use such 
measurements as input for a practical algorithm with 
which to predict field intensities and correlations away 
from the source. Initially we consider radiation into free 
space in Sec. IIVI bv studying a simple model source and a 
more realistic source obtained from a full field simulation. 
In Sec. El we describe an application to a problem with 
reflecting boundaries, which is a first step towards our 
ultimate goal of extending the method to propagation of 
CEs in more complex environments such as cavities and 
larger structures. 

We start from a planar source at ^ = 0, parametrized 
by coordinates x = (xi, • • • Xd) with d = 1 or 2 in general, 
and radiating into in the half-space z > 0. We aim to 
predict the CE 

rz{xB,XA) = {'lp{xB,z)'tp*{xA,z)) ( 1 ) 

for z > 0 under the assumption than it can be mea¬ 
sured (or otherwise modelled) near the source screen 
z = 0. Here, (.) denotes an ensemble average over differ¬ 
ent source field correlations such as a time or frequency- 
band average. Eurthermore, '^(x, z) denotes one of the 
tangential field components in the frequency domain. 


The results easily extend to cross-correlation between dif¬ 
ferent components. 

In the past, the focus has often been on predicting 
the propagation of probability density functions (PDE) 
of waves passing through time-domain random [HI or 
turbulent [23| media. In our approach, the propagation 
itself is treated deterministically, whereas the radiation 
from the source is characterised statistically. This can be 
done, for example, by measuring the spatial field along 
a surface close to the source and determining the source 
CE by averaging the signal over time. We thereby elim¬ 
inate statistical fluctuations carried by the wave fields 
by ensemble averaging physical observables over suitable 
parameters. 

We now present the CE propagation rule explicitly for 
single field components. Polarisation effects can also be 
accounted for by propagating the field-field correlation 
tensor, which can be derived from the dyadic free-space 
Green’s function 

The field in the region 2 : > 0 is naturally presented in 
terms of the partial Eourier transform. 


/ 


^{p,z)= I e-^'^P--i,ix,z)x, 


where '0(x, z) denotes a field component on the screen 
itself {z = 0) or to its right {z > 0) and k is the wave 
vector. The radiated fields can then be reconstructed us¬ 
ing the evolution of this partial field. This can be calcu¬ 
lated by using the dyadic second Green identity which, in 
a source-free region, becomes the dyadic version of Huy- 
gen’s principle [HI- Being a convolution integral, the par¬ 
tial Eourier transform of the surface integral transforms 
to an algebraic equation. Then, the boundary conditions 
given by the fields sampled in the near-held region of the 
source can be used to eliminate the magnetic held in such 
an equation. The result of this procedure, restricted to 
the electric held components parallel to the source plane, 
is the following inhomogeneous plane-wave solution 

(t> (p, z) = (p, 0) , (2) 


where 


T{p) 


- ipp if < 1 

W\p\ ~ f if bl^ > f- 


( 3 ) 


Eor the moment, we neglect waves incident from the right 
and thus only describe radiation from a strong, direc¬ 
tional source; including incoming waves at the interface 
can be introduced formally using the boundary integral 
equations according to the discussion in [HI- An exam¬ 
ple of this scenario, involving a planar rehector beyond 
the source, will be given in Sec. |Vl Here, p = (pi,... ,pd) 
takes the meaning of a momentum tangential to the d 
dimensional source plane. In the ray-dynamical limit, we 
may identify 


|p| = sina , 

T (p) = pz = cosa , 


( 4 ) 

( 5 ) 




3 


where the angle a describes the direction of the ray with 
respect to the local outward normal to the source. In 
this perspective, T {p) represents a generalized kinetic 
energy of the ray. The case |pp > 1 in ([3|) corresponds to 
evanescent propagation, which does not contribute to the 
far-held, but may be detectable in the near held; see also 
the discussion in Secs. lIVBi II V Cl In order to represent 
wave helds in phase-space using canonical coordinates 
(x,p) parallel to the source plane, we dehne the WDF 

Wz{x^p)= f s/2,x — s/2) s = 

fkVf ' 


the source, that is, Vq{xb^xa) = ro(x -h s/2^x — s/2) 
varies with x on a larger-than-wavelength scale. In that 
case, signihcant contributions to m are obtained only 
for small q and we can expand the phase diherence 
AT (p, q) = T{p q/2) — T'^{p — q/2) around q = 0. 

In the region |pp < 1 corresponding to propagating 
waves, the difference AT receives contributions only from 
odd powers of q. Neglecting cubic and higher order terms 
we hnd that 



This is the Frobenius-Perron (FP) propagator (SOlfor 
radiation into free space and leads to the evolution 


Upon insertion of (|2]) in m, and by exploiting the inverse 
transformation to represent the source correlation (at z = 
0) in terms of the source Wigner function IFo {x^p)^ we 
hnd 

w, (x,p) = I g, ix,p,x',p') m) {x',p') x' p'. (7) 

This provides us with a propagator of the Wigner func¬ 
tion taking the form 

gz{x,p,x',p') = 5{p-p') (8) 

^ J ^ik{x-x')-g+ikz(T(p+q/2)-T-{p-g/2)) 

where the ^-function represents translational invariance 
in X and the corresponding conservation of momentum. 
Eq. dHI provides a scheme to propagate wave densities in 
phase-space for arbitrary sources, no matter how complex 
or rapidly varying. The propagation of the correlation 
functions themselves can subsequently be retrieved by 
an inverse Fourier transform of ([7]). That is, 

T^{xb,xa) = (^■^') J e''^'’'PW^{x,p) p, ( 9 ) 

where x = {xa +xb)/ 2, and s = xb — xa- The intensiW 
Iz as function of the distance 2; can be retrieved using [3 

h{x) =T^{x,x) = Jwz{x,p)p. (10) 


Wz{x,p)^Wo(^x-^yp') (12) 

of the WDF in the region |pp < 1. This approximation 
is equivalent to identifying the propagation of the WDF 
with the propagation of phase space densities along rays 
according to the evolution 


X 


zp 

W) 


p = p'. 


(13) 


The paraxial approximation is obtained by using the lin¬ 
earized flow in the regime |pp ^ 1 [6[. Outside this 
regime, the full flow has asymptotes along |pp = 1, the 
transition to evanescent propagation. Note that inserting 
Eq. (IT^ in Eq. (j9j) and including only the contributions 
from the propagating region |pp < 1 leads to 

{xb,xa) « J e^’^P^Wo /x - P 

In the paraxial regime, that is for |pp ^ 1, and for quasi- 
homogeneous sources, Eq. m retrieves the well-known 
van Cittert-Zernike theorem (VCZT) or generalisation 
thereof [nisi. In the next section we will show how to 
extend the VCZT to the non-par axial regime and includ¬ 
ing evanescent waves. 

Expanding AT further in q leads to an improved prop¬ 
agator which is capable of mapping less homogeneous 
CE’s. Including cubic terms leads in the 2D case, for 
example, to the Airy form 


III. RAY TRACING APPROXIMATIONS 

Asymptotic approximation of the propagator (0 leads 
to a direct propagation method for the WDE in terms of 
rays p, [23-0, • We will give a derivation of this ray 

limit below and will also discuss more subtle wave effects 
such as evanescent decay into the near-field and higher 
order (in 1/k) wave corrections. 

The simplest ray-based approximation is obtained un¬ 
der the assumption that the CE is quasihomogeneous at 


Gz {X,p,x\p') ^5a(^X-x' - 1 /^') ( 15 ) 

where 

^a{xi) = aM{au)^ 

with a = 2k/{kzT"\p))^^^ and Ai denotes the Airy func¬ 
tion. Note that limi^i^oo ^a{xi) = 6(u), so the EP form is 
obtained in the limit of large wavenumber k as expected. 
Similar results have been obtained in the context of the 










4 


j^ro pagation of EM waves through inhomogeneous media 

Further improvements over the basic FP propagation 
m are obtained by accounting for evanescent decay into 
the near-held, which emerges from contributions |pp > 1 
in Eq. (0. Since the kinetic operators in Eq. (0 now add 
constructively, the leading contribution is formed by the 
zeroth order term in the expansion of AT(p, g) and we 
obtain 


Wz{x,p)^e 1 (x,p) , |pp > 1. (16) 

Improved approximations may be achieved by treating 
the exponent beyond leading order, but we hnd that (iroi) 
gives a good description of evanescent decay already as 
discussed in the next section. 



X [m] 

FIG. 1: Magnitude of the WDF of a ID Gauss-Schell corre¬ 
lation function. The radiation frequency is / = 1 GHz corre¬ 
sponding to A = 0.3 m. The distribution widths are cr^ = 1.0 
m, and as = 0.1 m. 


IV. RADIATION INTO FREE SPACE 


We now test the effectiveness of the FP propagator in 
the simple case of radiation into free space. The hrst ex¬ 
ample treated in Sec. II V Al assumes a quasi-homogenuous 
source distributed according to the Gauss-Schell model. 
We will examine the near-held behaviour in more de¬ 
tail in Sec. lIVBi considering in particular the limit of 
completely uncorrelated sources. In a second example in 
Sec.IVCl we consider a more complex set-up mimicking 
the more realistic sources expected in typical EM appli¬ 
cations. We will restrict ourselves in these examples to 
2D models (so d = 1) and characterise the behaviour of 
held-held correlations by focusing on the propagation of 
one held component along z. We select the held tangent 
to the source. 


A. Propagation of Gauss-Schell model 


Using a simple 2D model for the emission of partially 
coherent EM radiation, we assume a source correlation 
in terms of a truncated ID Gaussian-Schell model 


Fo {xb.xa) = /oexp 


{xb - XAf 


(xb + XAf 


exp 

8cr2 


xXi {xb) Xi (xa) 
(17) 


where I is the length of the source. Here, the character¬ 
istic functions 


»<*>={«: n 4 ’ 

account for the hnite size of the source. The quasi¬ 
homogeneity condition can be expressed through de¬ 
manding as ^ X ax, where A = 27r/k is the optical 


wavelength. The source WDF is then found to be 


Wo {x,p) = Jo exp 

erf ( 


k‘^p‘^a‘1 
-a, exp I-^ 


sV2 ^/2 


X 

2a‘i 


_ erf ( 


,y/2 y/2 


(19) 


For extended sources, for which / ^ A, and for x inside 
the region occupied by the source, Eq. m simplifies to 

Wo {x,p) ^ V^aglo exp 

Fig. [1] shows the WDF in phase-space at 2) = 0 for a 
spatially extended source [^. Here, and in all other 
computations with the Gauss-Schell model, we work at a 
frequency of operation of 1 GHz corresponding to A = 0.3 
m and choose da, = 1.0 m, as =0.1 m. 

Fig. [2] shows the propagation of ([T9]) as computed 
through the full integral operator ([7|) together with the 
propagation obtained by the Frobenius-Perron approx¬ 
imation m- There is surprisingly good agreement be¬ 
tween the exact and approximate behaviour even far from 
the paraxial regime. This is remarkable given that the ray 
tracing approximation is only valid to leading order. This 
constitutes a major computational advantage as the FP 
approximation reduces an integral equation to a simple 
coordinate transformation. The overall behaviour shown 
in Fig. |2ja) and Fig. |2jb) reflects the distribution shear¬ 
ing due to the geometrical ray propagation based on Eq. 
(HI; see also p. The CFs can now be obtained by a 
back transformation according to Eq. © and are shown 
in Fig. |2]^c) and Fig. |2fd). 


x^ 

2al 


k^p^ai 


( 20 ) 


B. Non-paraxial Van Cittert-Zernike theorem 

In the following, we will focus on near-held effects for 
small distances from the source as a function of the source 
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FIG. 2: Propagated Wigner distribution function (a), (b) and 
correlation function (c), (d) of a partially-coherent (ID Gauss- 
Schell) near-homogeneous source: exact ((a) and (c)) vs ap¬ 
proximate - Frobenius-Perron - ((b) and (d)) computation at 
z = 10A. The radiation frequency is / = 1 GHz. The distri¬ 
bution widths are (Jcc = 1 m, and cts = 0.1 m. 



FIG. 3: Evolution of the correlation function along the line 
X = 0 for a selected (partially) correlated source, from 
kz = 0.42 (blue solid line) to kz = 1.68 (brown solid line). 
In the presence of evanescent waves the correlation width in¬ 
creases rapidly before the correlation function becomes a sine 
function, whose width then increases linearly according to the 
VGZT. 


correlation parameter (jg • We are in particular interested 
in how the correlation length propagates in the near- field 
before reaching the linear VCZT regime. 

In the near-held limit, the WDF shows exponentially 
decaying evanescent components according to m, while 
the WDF remains essentially unchanged for the propa¬ 
gating part |pp < 1. This leads to a model for the WDF 
with source distribution (f2Q|) of the form 


Wz{x,p) ^ V^CTs/oexp 


2crJ 


k^p^af 


1 if |pp < 1 

e-2fc^Vlpl"-i if |p|2 > 1 _ 


( 21 ) 


Far enough from the source, such that evanescent compo¬ 
nents have completely decayed, while close enough that 
evolution in the propagating region of phase space can 
still be neglected, we model the WDF using 


A. The second moment of the correlation function is not 
dehned and cannot therefore be used to dehne a corre¬ 
lation length. Instead, we dehne the correlation lengths 
to be the spacing at which the correlation has fallen by 
a factor 1 / 

Tz{x + As/2, x-\- As/2)/Tz{x,x) = (23) 

Note that for a Gaussian correlation function such as 
assumed for the source in da, this dehnition coincides 
with the standard variance: As = ag. 

With the dehnition adopted above, we can now obtain 
the correlation lengths from exact wave propagation cal¬ 
culations. The results are shown in Fig. HJas a function 
of the distance z for diherent source correlation lengths 
ag. The universal regime is shown as the blue dashed 
line. The VCZT regime starts around kz > 1. 

From Eqs. m and (j22p . one can estimate the growth 
rate both in the near and the far held. Including non¬ 
par axial ehects, there are three diherent regimes: 


Wz{x,p) ^ v^cTs/oexp 


k/^p^cr/^ 


X 

2 

J 1 if |p|2 


0 


< 1 

if \p\^ > 1. 


( 22 ) 


Using the inverse Fourier transform, Eq. m, we now ob¬ 
tain the correlation function from the WDF given by Eq. 
([21]) or Eq. ([22|) acccordingly. In Fig. [3] we show the 
resulting near-held evolution of the correlation function, 
placing the midpoint x = {xa -\r xb)/2. = 0 at the centre 
of the source. 

One observes that in the near-held regime, the width 
As of the correlation function is smaller than A, but it 
increases rapidly towards A as z approaches and exceeds 


i) in the deep near held, with kz < \ and kug < 1, 
the correlation length increases linearly with a slope 
that is independent of the frequency as well as of 
Gg and 

ii) a no-growth regime with As = const. 

— For kGg < 1, one hnds kAs ^ 1 in the range 
1 < kz < kGx] 

— For kGg > 1, then As = Gg in the range 0 < 

2) < kGx- 

The latter regime has already been described by 
Cerbino for paraxial sources |lq : 
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near field in this regime. This is important for sources 
that show fluctuations on scales smaller than the wave¬ 
length, such as in the case of a fully uncorrelated source 
(Js = 0, which may serve as a model for thermal sources 


The plateau behaviour corresponding to regime (ii) 
arises when ^ is sufficiently large that (|2^ describes the 
WDF, while kas 1. The correlation function is then 
proportional to the inverse Fourier transform 


^z{s) ^ —smc{ks) 


of the function 


FIG. 4: Evolution of the correlation length As with z on 
a log-log scale at x = 0 and for a partially correlated ID 
source with kas = 0.002 (cyan solid line)... 6.3 (red solid 
line). Evanescent waves drive the rapid spatial (blue dashed 
line) increase preceding the VCZT regime (red dashed line). 
A plateau is observed between kz ^ 1 and the Rayleigh range 
discovered in Em. 


W,{x,p) 



bP < 1, 
bP > 1 


(26) 


of p. In this case the correlation length defined by Eq. 
(|2^ takes the form 


kAs ^ 1.6443 


(27) 


hi) the VCZT regime for large z with a linear growth 
of the correlation length according to 

As (X , (24) 

with a slope depending on the ratio of wavelength 
to source dimension £. 

We now motivate these three regimes in more detail, 
beginning with case (i), which corresponds to kag < 1, 
kz < 1. The WDF Wz[x,p) described by ([21]) then de¬ 
cays slowly along the p axis as \p\ increases beyond the 
propagating region |pp < 1. In the extreme nearfleld the 
correlation function is proportional to the inverse Fourier 
transform of the function 

W,{x,p) ^ 


of p, that is. 


r.W - 


2kz/'K 

{2kzY -h {ksY 


The correlation length defined by Eq. (|23|) then takes the 
form 


As 



1.6109 z. 


(25) 


independent of (jg- It should be noted that if the condi¬ 
tion k(7g < 1 is breached, then the Gaussian decay in p 
present in (|2^ becomes the dominant feature and instead 
a limiting plateau level 


As = ag 

occurs, see Fig. (H Note that in this case the plateau 
extends all the way to z = 0 and the linear regime of 
case (i) is not seen. 

Finally, regime (iii) applies once evolution of the phase 
space takes effect in the propagating region |pp < 1. 
Assuming the quasihomogeneous case ax ^ I and con¬ 
sidering first only the near-held region kag kz 
then for a given midpoint x the finite size of the source 
reduces the support in p of the Wigner function and (j26|) 
is replaced by 

fl ^ r < , 

Wz{x,p) ~ < ’ Yz‘^ + {x-ll2y ^Jz‘^ + {x+l/2y ’ 

[O, otherwise. 

For simplicity consider the case x = 0. Then the correla¬ 
tion function obtained from the inverse Fourier transform 
of this function is 

r^(s) - Isinc ( ^ 

^ (yi + (2V0V 


That is, we And in regime (i) that evanescent decay of 
the sub-wavelength correlations in the source dominates 
in such a way that there is a universal growth rate in the 
correlation length. The numerical value of the slope in 
([25]) is particular to the form taken in Eq. (|23]) for the 
correlation length, but the qualitative conclusion applies 
more generally. The presence of evanescent waves thus 
leads to a rapid increase of the correlation length in the 


and the correlation length defined by ([23]) takes the form 

As fa 0.2617Av'l + (2V0^ (28) 

generalising (ITT]) . In the farfleld z ^ we And 
As fa 0.2617 x ^ 
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FIG. 5: Set of parallel metallic wires running above an oblique 
perfect electric conductor (PEC) ground plane. This complex 
source emits radiation in the half-space z > 0. The full-wave 
TLM simulation has been carried out for the configuration 
hi = 0.3025 m, hu = 0.358 m, / = 1 m, d = 0.06 m. 

T 

(where the numerical prefactor is particular to the con¬ 
vention (ESI)). Alternatively, if ax ^ /, then the screen 
length becomes unimportant and ax provides the length 
scale appropriate to the source intensity. An analogous 
calculation then allows us instead to recover the basic 
form 


ax 

of the VCZT for ax- 


C. Application to a complex source 

The field at the source at 2 : = 0 is often produced by a 
complex process such as tracks on a printed circuit board 
or integrated circuits in electronic devices; the radiation 
produced in the source region then propagates into free 
space. We model such a complex source here by a set of 
N metallic wires driven by random time-domain voltages, 
as illustrated in Fig. 0 a realization of the voltage Sp (t) 
driving a pin of the bundle is reported in Fig.[6]along with 
its spectrum Sp (/). This represents a typical problem 
in EMC where one tries to obtain statistical information 
about an erratic signal. 

The presence of a perfect electric conductor along an 
oblique plane makes the source radiate only into the half¬ 
space z > 0: this mimics a configuration that is widely 
used in the design of printed circuit boards. We use N = 
12 wires very close to each other and to the metallic plane 
in terms of wavelength. Therefore, it is reasonable to 


think of this circuit as a collection of random sources of 
partially coherent radiation. 

The exact fields emitted from such a complex structure 
are computed through an in-house Transmission Line 
Matrix (TLM) code [^. This is a time domain method 
for modelling 3D electromagnetic field interactions with 
complex structures that may include a variety of materi¬ 
als. The technique is based on the equivalence between 
electric and magnetic fields and the voltages and currents 
on a network of transmission lines. After discretizing 
space, the fields in individual cells are modelled by trans¬ 
mission lines incident from each cell-face and intersecting 
at the cell centre forming a junction. Each of these or¬ 
thogonal transmission lines allows for the propagation 
of electromagnetic waves. The waves are characterised 
by voltage and current and their associated electric and 
magnetic fields. In order to obtain the desired correla¬ 
tion functions, we sample the numerically obtained fields 
in a plane above the tracks at different times in order to 
create a suitable ensemble of uncorrelated circuit reali¬ 
sations. This is used as a basis for calculating field-field 
correlation functions and their Wigner functions both in 
the near- and far-held. 

Figs. [ 7 ] (a) and (b) show the comparison between the 
WDF as computed through the full-wave (TLM) simula¬ 
tions, and the WDF obtained by the FP approximation 
([12]) in the far-held at 2 : = 2.3A. In the TLM calculation, 
the full time-dependent held is propagated out from the 
source, while in the FP approximation, the WDF ob¬ 
tained from the signal at the source (as shown in Fig. 
(6] see also Fig. [8] (a)) is propagated according to ([T^ . 
There is good agreement between the behavior predicted 
from full-wave simulations and the FP approximation, 
even though the source exhibits strong inhomogeneities. 
Interestingly, Fig. [ 7 ] shows the same Wigner distribu¬ 
tion shearing as in Fig. [2] following the geometrical in¬ 
terpretation m of the correlation propagation. It is 
worth stressing that such a Wigner function challenges 
the FP approximation ((Hj), whose underlying assump- 
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FIG. 9: TLM versus approximate energy density of the field 
radiated far from the source. 


FIG. 7: Far-field WDF (upper plots) and correlation func¬ 
tions (lower plots) at z = 2.3 A: comparison between TLM 
computation (left column) and Frobenius-Perron analytical 
approximation (|12l) (right column). The propagated correla¬ 
tion has been calculated through 



FIG. 8: (a) Source distribution; near-held WDF: comparison 
between TLM computation and evanescent-wave approxima¬ 
tion (ITbl) for (b) z = 0.1 A; (c) z = 0.2A and (d) z = 0.4A along 

X = 0. 

tion is quasi-homogeneity. Note that we can always also 
switch to the exact transport rule ©, which is computa¬ 
tionally more expensive than the FP approximation, but 
still orders of magnitudes faster than a full TLM calcu¬ 
lation. Propagated CFs as shown in Fig. 0 (lower plots 
(c) and (d)) are hnally obtained by applying the inverse 
Fourier transform 

Note that we also hud a pronounced broad side radia¬ 
tion around p ^ ±1 (corresponding to a ^ ±7r/2), and a 
strong asymmetry of the Wigner distribution due to the 
oblique metallic reflector. Those features can be captured 
by inspection of the WDF representation in phase-space, 
while they are less apparent in the propagated correlation 
function shown in Fig. 0(c) and (d). 

The source distribution Wo{x,p) as obtained from the 
radiated signal, see Fig. [6l is shown in Fig. [8] (a). Note 


that the region with |pp > 1 corresponds to evanescent 
contributions. In Fig. [8](b)-(d), a comparison between 
WDFs as computed through the full-wave (TLM) sim¬ 
ulations and those obtained using the WDF propaga¬ 
tor incorporating evanescent contributions, Eq. ([T6|) , are 
shown along the line x = 0. We find that propagation 
beyond z = 0.1 A results predominantly in an exponen¬ 
tial reduction of the WDF in the region |pp > 1. In the 
far-held, the radiation energy is restricted to the phase- 
space region \p\ < 1, as can be seen in the WDF in Fig. 
0(a) and (b). 

A comparison of (TLM) simulated and (Frobenius- 
Perron) approximate far-held propagated energy 

Ez = ^eo ^ eo J Wz {x,p) dp, (29) 

is shown in Fig. [9l We see that the two numerical meth¬ 
ods show qualitatively the same features, however, there 
are quantitative diherences. We think that these devia¬ 
tions are due to a diherence in the numerical treatment 
of the boundary conditions at x = ±3 m. While the FP 
approach has no difficulties in treating these boundaries 
as completely open, the TLM method needs to model 
this with absorbing boundary conditions. These condi¬ 
tions tend to be still slightly rehective, as is evident from 
the source distribution in Fig. [8] (a) around x = ±3, 
p = ±1. This comparison highlights another advantage 
of the Wigner function propagation method. The evalu¬ 
ation of the WDF of actual circuits can be done for the 
full electromagnetic field by using the approach described 
here component by component. 

V. REFLECTION OF PARTIALLY 
CORRELATED SOURCES 

Having developed a framework for the propagation of 
CFs in free space, we are now interested in tackling the 
case of reflection from planar boundaries. In particular 
we would like to test the FP approximation in the pres¬ 
ence of interference. It is then interesting to solve the 
canonical situation depicted in Fig. [TOl where a planar 
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Reflector 


FIG. 10: Arbitrary planar electromagnetic source emitting in 
the half-space z > 0, in presence of a planar metallic bound¬ 
ary. 

reflector is located at distance z — L from the source at 

z = 0. 

The reflecting boundary is here for simplicity assumed 
to be parallel to the source plane, indefinitely extended 
in the a%-plane, and made of an ideal perfect electric con¬ 
ductor (PEC). Therefore, for electric (TE) or magnetic 
(TM) fields perpendicular to i, the Eresnel reflection co¬ 
efficient reads r (a) = —1, for all incoming angles a [38^ . 
We again consider for simplicity only a scalar field, or 
a single component of the vector field, emitted from the 
source. 


A. Theory 

Consider a plane located at an arbitrary longitudinal 
coordinate z = D between source and detector. The 
field distribution in the plane consists then of two con¬ 
tributions: the direct wave coming from the source, and 
the reflected wave bouncing off the reflector back to the 
source, that is, 

^ (p, z) = Q) _ ^ikDT(p)+i2kAT(p)^^p^ ^ 

(30) 

where (j) (p, 0) is the field at the source plane 2 : = 0, T{p) is 
defined as in ([3]), and A = L — D. The momentum space 
CE is formed as the product of the two fields in dsni) and 
an ensemble average is taken as in Eq. (HJ. By plugging 
the closed-form expression (I30II into the definition of the 
WDF dSl), we find the phase space representation 

Wd ( p , x) = Wd ( p , x) + W2L-D (p, x) - [Wa ( p , x) -\r cc], 

(31) 

where the first two terms are direct and reflected con¬ 
tributions respectively, coming to the detector straight 
from the source or through the reflector, and the last 
two terms express the interference between direct and 
reflected waves with cc standing for the complex conju¬ 
gate. 

Eollowing the procedure described in the previous sub¬ 
section, it can be shown that direct and reflected terms 
in m can be calculated through the free-space propaga¬ 
tion scheme in ([7]) and (|8]), with z = D and z = 2L — D 
respectively, while the interference terms lead to 

Wa {p,x) = jj Qa {x,x'-,p,p') Wq {x',p') x'p' , (32) 


with a modified Green integral operator 

Qa{x,x';p,p') = d{p-p')-^ ( 33 ) 

(27r)^ 

^ J ^ik{x-x’)q+ikD{Tip+i)-T*ip-i))-i2kAT^p-i) ^ ^ 

Eor the class of statistically quasi-homogeneous sources, 
we may again expand the exponent in (|33|) in a Taylor 
series in q, and retain only terms up to first order. This 
results in a Erobenius-Perron approximation of the in¬ 
terference terms, leading to a phase-factor of the optical 
length A besides the Dirac’s delta in (ED. Adopting the 
same linear approximation for each term in (l3T]) gives the 
updated WDE 

WD{x,p)^Wo[x-^yp) ( 34 ) 

-2cos(2fcAT(p)) IFo {x - ^,p) . 

Similar expressions have been found in quantum mechan¬ 
ics [39[ and optics [5[ for two overlapping wave-functions. 

Again, the propagated CE can be obtained by the in¬ 
verse Eourier transform ([9j) of (|3T]) or the latter 

being closely related to the free-space VCZT. 

B. Numerical results 

We chose again an initial correlation density dis¬ 
tributed according to the Gauss-Schell model, Eq. (El, 
with corresponding source WDE shown in Eig. [H We 
work as usual at a frequency of operation of 1 GHz cor¬ 
responding to A = 0.3 m and choose ax = 1.0 m, ag = 0.1 
m. 

We further suppose a metallic mirror at L = 1.8 m 
(6A). The propagation of the intensity from the source 
to the mirror can be found by evolving the source WDE 
with the exact rule composed of Eqs. 0,0 and those 
for the interference terms, Eqs. (ED - and then 

inverse Eourier transforming the propagated WDE ac¬ 
cording to Eq. (|9]). The coherent energy Iz{x) reaching 
the scan plane at 2 : = D is given by Eq. m, that is, by 
considering the correlation function at s = 0. 

Eigure [TT] shows the behavior of the intensity Iz{x = 
0) near the mirror, from D = 1.0 m to D = 1.8 m. 
The solid black line is computed through the full Green’s 
integral operators (ED and (ED, while the dashed red 
line is obtained by the Erobenius-Perron approximation 
(l34|) . The oscillatory behaviour in (ED is due to the 
interference terms in the WDE. 

In Eigs.[T2]and[T3l we show the magnitude of the WDE 
and the associated CEs at a distance A = 1.25A (position 
A in Eig. fTT]) and at a distance A = 2A (position B in Eig. 
n]) from the mirror, respectively. While good agreement 
between the exact and the approximate propagation us¬ 
ing the EP approximation is achieved at position A, a 










10 



FIG. 11: Interference pattern formed by the intensity G 
(x = 0) along the longitudinal direction as the scan plane 
approaches the reflector. Exact (black solid line) versus ap¬ 
proximate (red dashed line) computations are compared. 


maximum in the correlation function, the same is not 
true at position B. Here the intensity is suppressed due 
to destructive interference and the magnitude of the CF is 
itself only of order 0{l/k). To obtain the good agreement 
shown in Fig. [131 we need to take into account higher or¬ 
der corrections in the WDF propagator such as using the 
Airy function integral kernel, Eq. (USD. The improvement 
when going from the leading order FP to the Airy func¬ 
tion approximation is shown in Fig. [13] (b) to (c), which 
need to be compared with the exact WF Fig. [13] (a); 
the corresponding propagated CF is displayed in Fig. [13] 
(d). Only after going beyond the FP approximation in 
this way are we able to reconstruct the fine structure of 
the WDF. This finding is not surprising, but remarkable 
nevertheless; computing WDFs in a multi-scattering en¬ 
vironment will encounter exactly these problems and we 
have shown that the Airy-function approximation - still 
faster than a full WDF propagation - can handle inter¬ 
ference corrections successfully. We note that these cor¬ 
rections have been reported also in the “diffusive” Green 
function presented in j^ . 


VI. CONCLUSION 



-10 0 10 


X [m] 



-10 0 10 
X [m] 


FIG. 12: Magnitude of the WDF of a ID Gauss-Schell source: 
exact (left plots) versus approximate (right plots) computa¬ 
tion at A = 1.25A (position A in Fig. fTTl) . Related correlation 
functions are reported in the lower plots. 



FIG. 13: Magnitude of the WDF of a ID Gauss-Schell source: 
exact (a) versus approximate using only FP approximation 
(b) and using the Airy approximation (c) at A = 2A (position 
B in Fig. 1111) . The correlation function corresponding to (c) 
is reported in (d). 


An exact propagator has been derived for field-field 
correlation functions of complex sources. It has been ap¬ 
plied to a problem mimicking EM radiation from a com¬ 
plex source; extending this to other wave problems such 
as in vibro-acoustics or quantum mechanics is straight¬ 
forward. The phase-space representation based on the 
Wigner function provides a useful means of physically in¬ 
terpreting the propagated data. It also serves as a very ef¬ 
ficient computational technique both for an exact propa¬ 
gation of CFs and in terms of a ray approximation leading 
to the Frobenius-Perron operator. This provides a good 
description of the propagated data even when applied 
to source data that are relatively far from homogeneity. 
Where necessary, more heterogeneous sources can be ac¬ 


counted for by higher-order approximations leading to an 
Airy propagator. This propagator proved important in 
the case of a planar random source emitting in presence 
of a planar reflector, for which we are able to reconstruct 
the fine structure of the phase space in presence of in¬ 
terference. Evanescent decay into the near field can also 
be accounted for using simple propagation rules. These 
rules has been used to investigate the effect of evanescent 
waves in near-held correlation functions. For source cor¬ 
relations exhibiting smaller-than-wavelength scales, we 
predicted a rapid initial increase of the correlation length 
(with distance from the source), before it saturates with 
the onset of the Van Cittert-Zernike behaviour at a dis¬ 
tance of a wavelength. The approximations used have 




















been validated through full-wave simulations using model 
sources and numerical sources exhibiting strong statisti¬ 
cal inhomogeneities. 
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